clear
clc

load data
X = Features;
Y = Label;
clear Features Label

%%

data_centered = bsxfun(@minus, X, mean(X, 1));
[~, score_centered] = princomp(data_centered, 'econ');
score_centered = score_centered(:, 1:2);
% Visualize data centered
figure(1)
subplot(121)
hold off
gscatter(score_centered(:, 1), score_centered(:, 2), Y)
title('Data centered PCA')
xlabel('PC1')
ylabel('PC2')
axis tight

%%

data_normed = zscore(X);
[~, score_normed] = princomp(data_normed, 'econ');
score_normed = score_normed(:, 1:2);
% Visualize data normalized
subplot(222)
hold off
gscatter(score_normed(:, 1), score_normed(:, 2), Y)
title('Data normalized PCA')
xlabel('PC1')
ylabel('PC2')
axis equal
axis tight

%%

aplsc_model = aplsc(data_normed, Y);
n_dim = 1; % Could be set to integer greater than 1
Y_hat = test_aplsc(aplsc_model, data_normed, n_dim);
[accuracy, sensitivity, specificity] = calc_classification_performance(Y, Y_hat);

%%

E = data_normed;
t1 = E*aplsc_model.W(:, 1);
E = E - t1*aplsc_model.P(:, 1)';
t2 = E*aplsc_model.W(:, 2);

subplot(224)
hold off
gscatter(t1, t2, Y)
title('Data normalized APLSC')
axis equal
axis tight

m1 = aplsc_model.V(1, 1)*aplsc_model.Q(1);
m2 = aplsc_model.V(2, 2)*aplsc_model.Q(2);
x = -5:0.02:5;
y = -m1*(x - aplsc_model.b)/m2;
hold on
plot(x, y, 'k') % The discriminant hyperplane
